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Dynamic Point-Formation in Dielectric Fluids 
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We use boundary-integral methods to compute the time-dependent deformation of a drop of 
dielectric fluid immersed in another dielectric fluid in a uniform electric field E. Steady state theory 
predicts, when the permittivity ratio, /3, is large enough, a conical interface can exist at two cone 
angles, with #<(/?) stable and 6>{/3) unstable. Our numerical evidence instead shows a dynamical 
process which produces a cone-formation and a transient finite-time singularity, when E and j3 are 
above their critical values. Based on a scaling analysis of the electric stress and the fluid motion, 
we are able to apply approximate boundary conditions to compute the evolution of the tip region. 
We find in our non-equilibrium case where the electric stress is substantially larger than the surface 
tension, the ratio of the electric stress to the surface tension in the newly-grown cone region can 
converge to a (5 dependent value, a c (/3) > 1, while the cone angle converges to #<(/?). This new 
dynamical solution is self-similar. 

PACS numbers: 47.11. +j, 47.20.-k, 68.05.-n 



The formation of conical ends on fluid-fluid inter- 
faces in strong electric/magnetic fields has been seen 
in various electrospraying and ferrofluid experiments 
0H II S Building on the work of Taylor @, Li 
et al. and Ramos & Castellanos studied the electrostat- 
ics of an infinite cone with semi- vertical angle Oq formed 
between two dielectric fluids with permittivity ratio (3. 
In spherical coordinates, the electric stress tr e ~ r 2 ^ -1 '. 
In an equilibrium cone, this stress must be balanced by 
the surface tension, so that v must be 1/2. According to 
their analysis there are two such solutions of 6 , 

0<(/3) and 0>((3), which will occur for (3 > (3 C = 17.59. 
The former is said to be stable, the latter unstable 0. 

In contrast, this letter describes a dynamical fixed 
point in which a cone is formed transiently. At the fixed 
point, the cone angle is 0<(/3) so that surface stress and 
electric stress have the same scaling in the cone, with the 
ratio of the two being constant, but different from unity. 
The electric stress is the larger of the two, so that the 
total surface stress always acts to elongate the pointed 
region. 

We compute the time-dependent deformation of a di- 
electric drop (fluid 1) freely suspended in another dielec- 
tric fluid (fluid 2) in a uniform electric field. Both fluids 
are incompressible and have the same viscosity rj. There 
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is surface tension with coefficient 7 between the two flu- 
ids. The drop is axially symmetric and has round tips, 
with its shape represented by the radius function h(z) in 
cylindrical coordinates (r, z). r* denotes the radius of 
curvature at the tip. An electric field with strength E 
is applied in the z direction (Fig. Suppose the ini- 
tial radius of the drop is a. Respectively we use a, 7/77, 
7/a, (7/ae £2) 5 and (7£ /ae2) 5 to scale length, velocity, 
stress, electric field and surface charge density [jj. 

Following Sherwood , we study a situation in which 
Reynolds number is small so that the fluid flows via 
Stokes equation and the charge distributions are deter- 
mined by electrostatics. The surface charge density p can 
be expressed in the form of a boundary integral equation 
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where p(x) is the surface charge density at (h(x), x), g is 
a Green function; [3 denotes the permittivity ratio £1/62; 
L y is the range of z axis occupied by the drop [j| . From 
the surface charge density p, the normal and tangential 
component of the electric field can be calculated to ob- 
tain the jump in the electric stress across the interface. 
Sherwood also uses a boundary integral to determine the 
interface velocity 
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FIG. 1: A drop of dielectric fluid freely suspended in another 
dielectric fluid in a uniform electric field. 



where i and j refer to the zorr component, fj(y) is the 
j component of the total surface stress, and G denotes 
a Green function [ljj. The velocity, u, is then used to 
update the interface position. In our simulation, we ap- 
ply a boundary element method with many details sim- 
ilar to that described by Sherwood. We distribute mesh 
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FIG. 2: Development of a finite singularity at j3 = 18.5 > j3 c 
and E 2 = 0.410 > E 2 . (A) The initial shape (the equilibrium 
shape at E 2 = 0.400), and the final shape calculated (r* = 
10~ 12 ). (B) Formation of conical ends. (C) Diverging velocity 
at the tip. 
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FIG. 3: The shape and stresses at r* — 10 -8 (dashed with 
dots) and r* = 10 -10 (solid) from the simulation in Fig. 
(A) Slopes on the interface. (B) Ratios of the electric stress 
to the surface tension on the interface. Tan 8 > — 0.689654 
and Tan6 < = 0.469704 at /3 = 18.5. The later curve matches 
the earlier one on the right hand parts of the plots. 



points in proportion to the local curvature, and use a 
cubic spline to interpolate the interface between mesh 
points, a quartic polynomial to interpolate the surface 
charge density. The derived linear algebraic equations 
are solved by LU decomposition. A fourth-order Runge- 
Kutta scheme is applied to update the interface position. 

There exists a critical electric field E c ((3) for (3 > (3 C . 
When E < E c , the drop can reach equilibrium with 
round tips. We start our simulation from a sphere and 
apply a sufficiently small electric field. If the maximum 
velocity on the interface decreases to a value below 10~ 4 
following an exponential decay, we consider that the drop 
will reach equilibrium. After equilibrium is reached by a 
numerical extrapolation, the field is increased by a small 
amount. Through increasing the electric field step by 
step we find the critical electric field. 

When E > E c , a finite time singularity develops. From 
now on, we use (3 = 18.5 as our example, which has E 2 = 
0.4085 ± 0.0003. For instance, when we choose the initial 
shape to be the equilibrium shape at E 2 = 0.400 < E 2 



and suddenly apply E 2 — 0.410 > E 2 , the drop forms 
conical-like ends. The velocity at the tip dramatically 
increases as a critical time is approached (Fig. |2J. Here 
we can at most obtain about twelve decades of data in 
one calculation, due to the increasing number of mesh 
points required and roundoff error. 

Figure |31 shows how the shape and stresses evolve as 
the finite time singularity develops. The slope plot sug- 
gests we can partition the interface into three regions, the 
tip region, the conical-like region and the macroscopic re- 
gion. The conical-like region is the intermediate region 
with a small variation in the slope. Figure CHA) shows 
as r* decreases, the macroscopic region and the estab- 
lished part of the cone region almost remain intact, while 
part of the interface which used to be in the tip region 
now grows conical-like. This shows that in the course of 
r* — > 0, only the tip region changes in time, while the es- 
tablished part of the cone region nearly remains indepen- 
dent of time. A careful examination of Fig. [21 A) shows 
the conical-like region is not exactly a cone, because the 
slope of the newly-grown cone changes as r* decreases. 
As we shall see in more detail later, this slope approaches 
the value set by 8 < . Figure E^B) shows the ratio of the 
electric stress to the surface tension is larger than one 
in the tip region and the conical-like region. The stress 
ratio also does not change in the part of interface whose 
shape remains as the conical-like region expands. Hence 
the numerical evidence says that the shape of each part 
of the almost conical region and the stress ratio within 
that part remain frozen as the tip gets smaller. How- 
ever, as r* changes the slope and the stress ratio of the 
newly-grown part change too. So Figure and [3] may 
show an approach to a fixed point, but they do not show 
a fixed behavior themselves. Because there is a slow and 
not-quite uniform convergence to a fixed point, it is hard 
to estimate the critical time, t c , from our raw data. For 
this reason, we shall henceforth plot our results against 
tip radius, r* instead of trying to use t c — t. 

A scaling study, sometimes called an order of magni- 
tude analysis, enables us to estimate the sizes of different 
contributions, it indicates how the different regions affect 
one another. 

These estimates show that except for a uniform advec- 
tion, the stresses in the tip region determines the flow 
within that region and the subsequent shape of the tip. 
Specifically, the deformation of the tip region is domi- 
nantly caused by the local stress jump. The axial strain 
rate defined as du z (x)/dx measures how fast the inter- 
face deforms due to the axial velocity. Using , we can 
express the contribution to du z (x)/dx from the three re- 
gions. Respectively the tip region, the conical region and 
the macroscopic region have length scales r* , h(z) and 
1 (r* <C h(z) <C 1). The electric stresses in the three 
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regions have an order of magnitude 
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with < v < 1. A similar result applies to the surface 
tension a s , but with v = 1/2. An argument like that of 
Lister and Stone shows that the forces in the intermedi- 
ate and macroscopic region simply advect the tip region 
without significant contribution to the strain rate 

A followup study |l2| will describe in more detail how 
the scaling analysis of the electric stress works. For the 
present purposes, it suffices to say that the shape in the 
tip region mostly determines the electric stress within 
that region, except for a coefficient which only depends 
on the shape in the other regions and the applied electric 
field. If we change the shape in the other regions, the 
electric stress in the entire tip region will be changed by a 
factor which is independent of the shape in the tip region. 
Changing the applied electric field will have the same 
effect. So after we reshape the rest part of the interface, 
we can restore the electric stress in the tip region by 
applying a different electric field of certain strength. 

The scaling study permits us to construct approximate 
boundary conditions which then permits the accurate de- 
termination of the subsequent behavior of the tip. Basi- 
cally whenever we are about to run out of mesh points, 
we cut off the part of interface far away from the tip and 
replace it by a new shape profile which takes fewer mesh 
points. Then we restore the electric stress in the entire 
tip region, which we can accomplish by adjusting the ap- 
plied electric field to restore the electric stress at the tip 
to its value prior to the truncation. The tip regions of 
the prescribed new drop and the original drop will subse- 
quently evolve in the same way, because the deformation 
of the tip region is primarily driven by the local stress 
jump. We define the rescaled axial distance £ and radius 
function H(^) as 



€=(z- z t ip)/r* 



J7(0 - h(z)/r 



(4) 



We at least keep the part of interface with £ < I0 4 and 
typically match a spherical band to the center region, 
requiring the slope to be continuous at the truncation 
points. The center of the spherical band, which locates 
on the z axis, coincides with the center of the prescribed 
new drop. The error will be smaller if the truncation 
point is farther away from the tip. The "truncate and 
prescribe" idea was invented by Zhang and Lister |l3l |. 

Using the same initial condition as in Fig. |3 we calcu- 
late the evolution of the tip region for eighty decades of r* 
with approximate boundary conditions. We truncate the 
drop for twenty-four times, starting at r* = ICP 8 . Figure 
0] shows that the approximate boundary condition pro- 
duces the same result as the exact boundary condition 
without truncation at 10~ 12 < r* < 10~ 8 . The point 
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FIG. 4: Overlap between the results calculated with the ex- 
act boundary condition (solid curves) and the approximate 
boundary condition (dotted curves). The initial condition is 
given in Fig. [5] Evolution of the tip region is calculated for 
totally 80 decades of r*, only part of which are shown here. 
Out of the 80 decades of data, the curves at small r* con- 
verges as power laws in r* . The fit shows as r* — > 0: (A) 
The slope at £ = 100 converges to 0.47260 ± 0.00003, (B) the 
stress ratio at £ = 100 converges to 1.1302 ± 0.0004, and (C) 
dr*/dt converges to -0.01977 ± 0.00003. 
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FIG. 5: Convergence to a fixed point at /3 = 18.5. The 
initial shape and the applied field of simulation (i) in solid 
lines: a sphere, E 2 — 0.4748; simulation (ii) in dotted lines: 
a sphere, E 2 = 0.473; simulation (iii) in dashed lines: the 
equilibrium shape at E 2 = 0.408, currently E 2 = 0.478. 
The overlap between the approximate boundary condition re- 
sults and the exact boundary condition results is similarly 
checked like in Fig. H As r* -» 0, (A) H'(£ = 100) -> 
0.472595 ± 0.000005; (B) «(£ = 100) -> 1.13040 ± 0.00004; 
(C) dr*/dt -» -0.019780 ± 0.000003. 



with £ = 100 is pretty close to the cone region, so the 
slope and the stress ratio there can respectively reflect 
the angle and the stress ratio of the newly-grown cone. 
Out of the eighty decades of data, the curves at small r* 
can be adequately fitted as c + b(r*) p with p > 0. Thus 
H (£ = 100), a(£ = 100) and dr* /dt each approach lim- 
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FIG. 6: The self-similar solution at the fixed point a(£ = 
100) c = 1.130 at p = 18.5. The first 12 decades of data 
computed with the exact boundary condition in simulation 
(i) reveal: (A) H(£) is a constant for small r* , for example 
at £ = 10. (B) The intermediate region has a constant slope 
0.46970 nearly equal to Tan 0<(18.5) = 0.469704. (C) The 
velocity at the tip increases logarithmically. (D) The stress 
ratio in the intermediate region is a constant 1.130 > 1. In 
(B)(D), dashed curves: r* = 10 -8 , solid curves: r* = 10 -12 . 



iting values as r* — ► 0. Just below we shall show those 
limits are the same for different initial conditions. 

Further simulations show there exists a fixed point be- 
havior: the stress ratio of the newly-grown cone con- 
verges to a fixed value larger than unity. At (3 = 18.5, 
a(£ = 100) will converge to 1.130 as r* — > 0, if it is 
close to 1 when the drop starts to develop conical ends, 
regardless of the initial shapes. For example in Fig. 
simulation (ii) in solid lines and simulation (iii) in 
dashed lines have different initial conditions. The dot- 
ted lines and dashed lines can all be excellently fitted by 
c + b(r*) p . As r* — > 0, the two simulations give the same 
limits = 100) c = 0.4726, a{£ = 100) c = 1.130 and 

(dr* /dt) c = —0.0198. As you may notice, we have ob- 
tained the same limits in the simulation in Fig. The 
fitted values of the power are very close to each other in 
all the simulations, and we get p = 0.013 ± 0.001. In 
simulation (i), we purposely choose a particular initial 
condition to let a(£ = 100) equal to cv(£ = 100) c at an 
early stage in the cone formation. We see that a(£ = 100) 
stably stays at a(£ = 100) c for many decades of r* with 
H (£ = 100) and dr* /dt equal to the limits obtained 
from the curve fitting. We have concrete numerical ev- 
idence that a(£ = 100) c = 1.130 is a stable fixed point 
at j3 — 18.5. This fixed point is finally approached in a 
power law of r* . And simulation (i) gives the solution at 
this fixed point. 

At this fixed point, the tip region is self-similar and the 
intermediate region is a cone with the cone angle 9 < (f3) . 
The first twelve decades of data calculated with the exact 
boundary condition in simulation (i) reveal the solution 



at this fixed point. We find the following properties: (a) 
The shape profiles of the tip region are self-similar af- 
ter we rescale them by r*, for example H(£) at £ = 10 
is a constant for small r* [Fig. EI A)], (b) The ever ex- 
panding intermediate region has a constant slope equal 
to Tan 0<(/3) [Fig. MB)}, (c) The velocity at the tip in- 
creases logarithmically [Fig. HJC)]. (d) The stress ratio in 
the intermediate region is a constant substantially larger 
than one, which we call it a c [Fig. HjD)]. The values 
of a c and a(£ = 100) c are very close to each other, (e) 
dr* /dt is a constant [Fig. EfC)], which indicates that r* 
scales like t c — t. This self-similar solution has some qual- 
itative similarity with a scaling solution found by Lister 
and Stone (lll |. 

We find similar fixed points at other values of (3 > 
[3 C such as 19.0. To summarize, our numerical evidence 
shows a cone with the smaller cone angle 9 < (f3) can be 
formed transiently in a non-equilibrium case where the 
electric stress is not balanced by the surface tension. The 
angle is approached as the stress ratio in the newly-grown 
cone region converges to a fixed value a c ((3) > 1. The 
dynamical solution at this fixed point is self-similar. 
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